Brownian motion and Diffusion#

import numpy as np
import scipy 

from numpy.random import normal, choice, uniform

import ipywidgets as widgets
import matplotlib.pyplot as plt
plt.style.use('fast')

%matplotlib inline

Mean square displacement (MSD) of random walker#

  • After time t number of steps (time) how far has random walker moved from the origin?

  • We quantify this by computing means square displacement. Where mean is computed over N number of simulated (observed) trajectories.

\[ MSD(n)= \Big\langle \big [ Z_n (t) - Z_n (0) \big]^2 \Big \rangle \sim t^{1/2} \]
  • Einstein developed a theory of diffusion based on random walk ideas and obtained a key equation relating mean square displacement to time \((t = n \Delta t)\)

\[MSD(t) = 2d D \cdot t^{1/2}\]

MSD of a 1D random walk#

def rw_1d(T, N):
    '''
    L: trajectory length
    N: Number of trajecotries
    returns np.array with shape (T, N) 
    '''
    
    # Create random walks 
    r  = choice([-1,1], size=(T, N))
    
    #Accumulate position
    rw = r.cumsum(axis=0)

    #Set initial position 
    rw[0,:]=0 
    
    return rw
T, N = 2000, 1000
rw = rw_1d(T, N)

t = np.arange(T)

R2 = (rw[:, :] - rw[0, :])**2 # Notice we subtract initial time

msd =  np.mean(R2, axis=1)    # Notice we average over N

plt.loglog(t, np.sqrt(msd), lw=3) 

plt.loglog(t, np.sqrt(t), '--')

plt.title('Compute mean square deviation of 1D random walker',fontsize=15)
plt.xlabel('Number of steps, n',fontsize=15)
plt.ylabel(r'$MSD(n)$',fontsize=15);
../_images/03_Diffusion_5_0.png

MSD of a 2D random walk#

def rw_2d(T, N):
    '''2d random walk function:
    T: trajectory length
    N: Number of trajecotry
    returns np.array with shape (T, N)
    '''
    verteces = np.array([(1,  0),
                         (0,  1),
                         (-1, 0),
                         (0, -1)])
    
    rw       = verteces[choice([0,1,2,3], size=(T, N))]
    
    rw[0, :, :] = 0
    
    return rw.cumsum(axis=0)
traj = rw_2d(T=10000, N=1000)
#Simulate 2D random walk
T, N = 10000, 1000
traj = rw_2d(T, N)

#Compute RSD 
dx = (traj[:, :, 0]- traj[0, :, 0]) 
dy = (traj[:, :, 1]- traj[0, :, 1]) 

R2     = np.mean(dx**2 + dy**2, axis = 1)   # notice how we averaging over N

fig, ax  = plt.subplots(nrows=2, figsize=(10,10))

t = np.arange(T) # time axis
ax[1].loglog(t, np.sqrt(R2), lw=3, alpha=0.5);
ax[1].loglog(t, np.sqrt(t), '--');

ax[0].set_title('2D random walker',fontsize=15)
ax[0].plot(traj[:3000, :5, 0], traj[:3000, :5, 1]);
ax[0].set_xlabel('X')
ax[1].set_xlabel('Y')

ax[1].set_xlabel('Number of steps, n',fontsize=15)
ax[1].set_ylabel(r'$MSD(n)$',fontsize=15);
../_images/03_Diffusion_9_0.png

Brownian motion#

  • The Brownian motion describes the movement of a particle suspended in a fluid resulting from random collisions with the quick molecules in the fluid (diffusion).

  • A small particle undergoes large number of molecular collisions when going from one step to another or after \(dt\) time. As a result each displacement over time \(dt\) can be viewed as a sum of ranomd collisions which can be approximated by a normal distribution via Central Limit Theorem.

  • Thus to simulate brownian motion we draw ranom displacements from normally distribution.

\[x(t+dt)-x(t)=N(0,\sqrt{2D dt})\]

We assume we have started at position \(\mu=0\) and our variance is given by \(\sigma^2=2Dt\) Where D is the diffusion coefficnets which is related to parameters of discree random walk as shown in the lecture.

\[x(t+dt)=x(t)+\sqrt{2D dt} \cdot N(0,1)\]

In the last step we re-wrote brownian motion equation in a convenient way by shifting normally distributed radnom variable by \(\mu\) and scaling it by \(\sigma\)

\[N(\mu, \sigma^2) = \mu + \sigma N(0,1) \]
def brown(T, N, dt=1, D=1):
    
    """
    Creates 3D brownian path given:
    time T 
    N=1 trajecotires
    dt=1 timestep
    D=1 diffusion coeff
    returns np.array with shape (N, T, 3)
    """
    
    n = int(T/dt) # how many points to sample
    
    dR = np.sqrt(2*D*dt) * np.random.randn(N, n, 3) # 3D position of brownian particle
    
    R = np.cumsum(dR, axis=1) # accumulated 3D position of brownian particle
    
    return R
  • Below we proceed to simulate continuous random walk in 1D-3D

  • We will plot trajecotires and distributions of brownian particle using interactive plotting via ipywidgets and holoviews/plotly interface.

R=brown(T=3000, N=1000)
print(R.shape)
(1000, 3000, 3)
@widgets.interact(t=(10,3000-1))
def brownian_plot(t=10):
    
    fig, ax = plt.subplots(ncols=2)
    
    ax[0].plot(R[:5, :t, 0].T, R[:5, :t, 1].T);
    
    ax[1].hist(R[:, 10, 0], density=True, color='red')
    ax[1].hist(R[:, t, 0], density=True)
    
    ax[1].set_ylim([0,0.1])
    
    ax[0].set_ylim([-200, 200])
    ax[0].set_xlim([-200, 200])
    
    fig.tight_layout()
import holoviews as hv
hv.extension('plotly')

plots = []
for i in range(10):
    
    plot = hv.Path3D(R[i,:,:], label='3D random walk').opts(width=600, height=600, line_width=5)
    plots.append(plot)
    
hv.Overlay(plots) 
rw_curve = [hv.Curve((R[i,:,0], R[i,:,1])) for i in range(10)]

xdist = hv.Distribution(R[:,10,0], ['X'], ['P(X)'])
ydist = hv.Distribution(R[:,10,1], ['Y'], ['P(Y)'])

hv.Overlay(rw_curve) << ydist << xdist

Diffusion Equation#

The movement of individual random walkers \(\leftrightarrow\) density of walkers \(\rho(\vec{r},t)\)

Diffusion equation:

Formulated empirically as Fick’s laws

\[\frac{\partial\rho}{\partial t} = \mathcal{D}\nabla^2\rho\]
  • This is a 2nd order PDE! Unlike equations of motion diff eq shows irreersibile behaviour

  • This one exactly solvable. But in general reaction-diffusion PDEs difficult to solve analytically.

  • Can solve numerically by writing derivatives as finite differences!

  • Can also simulate via random walk!

  • Diffusion coefficient \(D\), Units \([L^2]/[T]\)

Important special case solution (here written in 1d):

\[\rho(x,t) = \frac{1}{\sqrt{2\pi \sigma(t)^2}}\exp\left(-\frac{x^2}{2\sigma(t)^2}\right),\]

where \(\sigma(t)=\sqrt{2{D}t}\)

  • density remains Gaussian

  • Gaussian becomes wider with time

  • check that this is indeed a solution by plugging into the diffusion equation!

def sigma(t, D = 1):
    return np.sqrt(2*D*t)

def gaussian(x, t):
    return  1/np.sqrt(2*np.pi*sigma(t)**2) * np.exp(-x**2/(2*sigma(t)**2)) #
@widgets.interact(t=(1,100,1))
def diffusion(t=0.001):
    
    R = brown(T=101, N=1000)
    x = np.linspace(-20, 20, 100)
    
    plt.plot(x, gaussian(x, 1), '--', color='orange', label='t=0')
    
    plt.plot(x, gaussian(x, t), lw=3, color='green', label=f't={t}')
    
    plt.hist(R[:,t,0], density=True, alpha=0.6, label='simulation hist')
    
    plt.legend()
    plt.ylabel('$p(x)$')
    plt.xlabel('$x$')
    plt.xlim([-25, 25])

References#

The mighty little books

More in depth

On the applied side